R Notebook

Differential Expressed Genes

@see: http://dputhier.github.io/jgb71e-polytech-bioinfo-app/practical/rna-seq_R/rnaseq_diff_Snf2.html

Mediante las funciones del paquete de bioconductor “DESeq2”, se va a proceder a encontrar los genes que se expresan de forma diferente entre las muestras de los pacientes con glioblastoma de grado II y IV.

Paquetes necesarios

# Instalar, si no está instalado, DESeq2: herramienta para el análisis transcriptómico aguas abajo
if (!require("DESeq2", quietly = TRUE)) {
  BiocManager::install("DESeq2")
}

## Warning: package 'DESeq2' was built under R version 4.2.2

## Warning: package 'S4Vectors' was built under R version 4.2.1

## Warning: package 'BiocGenerics' was built under R version 4.2.1

## 
## Attaching package: 'BiocGenerics'

## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs

## The following objects are masked from 'package:base':
## 
##     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
##     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
##     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
##     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
##     Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
##     table, tapply, union, unique, unsplit, which.max, which.min

## 
## Attaching package: 'S4Vectors'

## The following objects are masked from 'package:base':
## 
##     expand.grid, I, unname

## Warning: package 'IRanges' was built under R version 4.2.1

## 
## Attaching package: 'IRanges'

## The following object is masked from 'package:grDevices':
## 
##     windows

## Warning: package 'GenomeInfoDb' was built under R version 4.2.2

## Warning: package 'SummarizedExperiment' was built under R version 4.2.1

## Warning: package 'MatrixGenerics' was built under R version 4.2.1

## Warning: package 'matrixStats' was built under R version 4.2.2

## 
## Attaching package: 'MatrixGenerics'

## The following objects are masked from 'package:matrixStats':
## 
##     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
##     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
##     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
##     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
##     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
##     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
##     colWeightedMeans, colWeightedMedians, colWeightedSds,
##     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
##     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
##     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
##     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
##     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
##     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
##     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
##     rowWeightedSds, rowWeightedVars

## Warning: package 'Biobase' was built under R version 4.2.1

## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.

## 
## Attaching package: 'Biobase'

## The following object is masked from 'package:MatrixGenerics':
## 
##     rowMedians

## The following objects are masked from 'package:matrixStats':
## 
##     anyMissing, rowMedians

Adquirir datos

Se adquieren los datos de la cunatificación generados con SALMON y importados con tximport. Por otro lado, también se adquieren los datos clínicos de los pacientes.

Previamente, se elimina los datos del pacientes 88, ya que por un error en el procedimiento del laboratorio, no es adecuado incluirlo en el análisis.

# --- Datos de expresion genética ---

# El objeto que contiene toda la información esta guardado como .RDS
txi <- readRDS("../import quant/txi.RDS")
# Trabajaremos con los counts
countData <- txi$counts

# No vamos a trabajar más con txi, así que como ocupa 35MB, lo elimino de la memoria
rm(txi)

# La matriz tiene, en los nombres de las muetras, un sufijo después de una barra baja, que es el "sample ID" que se utilizó para secuenciarlo en illumina. Para que no de problemas, se van a renombrar

# Elimino el sufijo "_Sxx" de los nombres de las columnas
lNewColnames <- unlist(strsplit(colnames(countData), "_"))[seq(1, 2 * ncol(countData), 2)]

# Asigno los nuevos nombres de las columnas a la matriz
colnames(countData) <- lNewColnames

# Muestro los primeros datos
# head(countData)

rm(lNewColnames)

# --- Datos clínicos ---

# Cargar la librería necesaria para leer archivos excel
library("readxl")

## Warning: package 'readxl' was built under R version 4.2.1

# Leer los datos del excel
clinicalData <- read_excel("../datos clinicos/dfClinical.xlsx")

# Muestro los primeros datos
# head(clinicalData)

# Los datos de la matriz estan ordenados como si los números fuesen carácteres, por lo que voy a transformar la variable "Sample" del dataframe de los datos clínicos a carácter, y después lo ordenaré de menor a mayor
clinicalData$Samples <- as.character(clinicalData$Samples)
clinicalData <- clinicalData[order(clinicalData$Samples),]

# Elimino la columna del paciente 88 de la matriz "countData" 
countData <- countData[,colnames(countData) != "88"]

# Elimino la fila del paciente 88 en los datos clinicos
clinicalData <- clinicalData[clinicalData$Samples != "88",]

# Nos aseguramos de que hay el número correcto de muestras
ncol(countData)

## [1] 27

nrow(clinicalData)

## [1] 27

# Mostramos
head(countData)

##                        103     105     112     113     131     151     165
## ENSG00000000003.15 280.690 269.872 437.320 354.990 565.889 474.571 393.377
## ENSG00000000005.6    0.000   1.000   2.000   5.000  14.000   0.000   1.000
## ENSG00000000419.14 468.499 657.489 417.199 670.997 476.536 303.273 464.651
## ENSG00000000457.14 292.789 289.795 250.316 338.842 339.367 258.275 297.864
## ENSG00000000460.17 148.557 101.415 201.240 132.638 132.400 169.704 270.774
## ENSG00000000938.13  37.000  84.008  82.001  36.000  87.000 108.000  39.180
##                         17     171     177     188     201     203     210
## ENSG00000000003.15 537.758 248.298 295.417 678.890 518.671 377.485 602.191
## ENSG00000000005.6    4.000   1.000   0.000   0.000   0.000   0.000   1.000
## ENSG00000000419.14 353.991 540.465 459.604 386.604 494.334 247.871 623.557
## ENSG00000000457.14 259.594 223.485 276.667 319.040 335.244 197.801 352.266
## ENSG00000000460.17 159.832 141.252 192.590 143.592 247.313 121.909 272.588
## ENSG00000000938.13  82.001  37.000  49.000  38.000 103.999 114.000  57.000
##                        218      30      37       5      63      64      70
## ENSG00000000003.15 636.489 418.424 410.440 597.728 460.811 463.123 625.235
## ENSG00000000005.6    2.000   3.000   0.000   0.000   1.000   1.000   1.000
## ENSG00000000419.14 616.128 881.789 460.925 457.024 566.649 324.236 852.521
## ENSG00000000457.14 370.288 464.448 238.741 257.118 287.889 256.463 359.007
## ENSG00000000460.17 302.957 212.938 222.348 178.380 145.989 125.019 226.644
## ENSG00000000938.13  18.000  88.000  70.999 144.001  44.000  59.000 144.590
##                         77      78      80      83      95      98
## ENSG00000000003.15 438.495 425.433 250.667 416.714 661.838 572.685
## ENSG00000000005.6    0.000   2.000   4.000   0.000   0.000   0.000
## ENSG00000000419.14 605.586 323.024 824.963 543.388 825.350 526.198
## ENSG00000000457.14 313.623 232.061 315.194 242.842 286.643 397.920
## ENSG00000000460.17 207.402 214.585 109.392 204.298 161.989 267.191
## ENSG00000000938.13  69.000  91.030 244.000  47.001  92.343  77.000

head(clinicalData)

Diseño del análisis

Mediante una función que facilita el paquete “DESeq2”, la información de ambos dataframes (abundancia de transcritos y datos clínicos) se “sumariza” en una matriz, donde las filas son los transcritos y las columnas los pacientes. En esta función, se debe de índican en base a que fáctor se quiere realizar el estudio. En este caso, se índica que se quiere realizar el experimento en base al fáctor “Grade” (el grado del glioblastoma).

# Cargar la librería
library("DESeq2")

# Aplicar la transformación que une counts + metadados
dds <- DESeqDataSetFromMatrix(countData = round(countData), colData = clinicalData, design = ~Grade)

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

# Mostramos la información del objeto (es un objeto especial de DESeq2)
dds

## class: DESeqDataSet 
## dim: 38244 27 
## metadata(1): version
## assays(1): counts
## rownames(38244): ENSG00000000003.15 ENSG00000000005.6 ...
##   ENSG00000291240.1 ENSG00000291292.1
## rowData names(0):
## colnames(27): 103 105 ... 95 98
## colData names(21): Samples Age.at.diagnosis ... fTumorLobule fSurvival

# Aplicamos un pre-filtrado, para eliminar las filas que tienen muy pocas lecturas.
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]
rm(keep)

# Guardamos el objeto resultante, para hacer uso de él en futuros análisis
saveRDS(dds, "./dds.RDS")

Diferential Expressed Genes

Aplicar la fución “DESeq” del paquete DESeq2 para obtener los datos del FC y el p valor. Estos datos permitirán aplicar un filtro para obtener los differential expressed genes o “DEGs” más relevantes.

library(DESeq2)

# Aplicar función "DESeq"
dds <- DESeq(dds, betaPrior=FALSE)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 774 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

resultsNames(dds)

## [1] "Intercept"      "Grade_IV_vs_II"

res <- results(dds)

head(res)

## log2 fold change (MLE): Grade IV vs II 
## Wald test p-value: Grade IV vs II 
## DataFrame with 6 rows and 6 columns
##                     baseMean log2FoldChange     lfcSE      stat     pvalue
##                    <numeric>      <numeric> <numeric> <numeric>  <numeric>
## ENSG00000000003.15 453.13289      0.0562706  0.279516  0.201314 0.84045267
## ENSG00000000005.6    1.56593     -0.0812530  1.530361 -0.053094 0.95765703
## ENSG00000000419.14 517.44826      0.4443501  0.252030  1.763081 0.07788683
## ENSG00000000457.14 291.63201      0.1676021  0.153023  1.095270 0.27339828
## ENSG00000000460.17 184.66633      0.7334334  0.281538  2.605092 0.00918496
## ENSG00000000938.13  77.75056      1.0953891  0.450076  2.433785 0.01494188
##                         padj
##                    <numeric>
## ENSG00000000003.15 0.9226214
## ENSG00000000005.6         NA
## ENSG00000000419.14 0.2498018
## ENSG00000000457.14 0.5024000
## ENSG00000000460.17 0.0783873
## ENSG00000000938.13 0.1022164

saveRDS(res, "./differentialExpresionResults.RDS")
write.csv(res, "./differentialExpresionResults.csv")

Doble filtro

Aplicamos el doble filtro (FC y p-valor) para obtener los genes más relevantes que se expresan de forma diferente.

# Omitimos los registros con NAs
resNotNans <- na.omit(res)

# Parámetros
# log2(FC) > 100
# p-value > 5
log2FC <- 2
pValue <- 0.05

# Filtro
filteredResults <- resNotNans[abs(resNotNans$log2FoldChange) > log2FC, ]
filteredResults <- filteredResults[filteredResults$pvalue < pValue, ]

filteredResults

## log2 fold change (MLE): Grade IV vs II 
## Wald test p-value: Grade IV vs II 
## DataFrame with 1302 rows and 6 columns
##                     baseMean log2FoldChange     lfcSE      stat      pvalue
##                    <numeric>      <numeric> <numeric> <numeric>   <numeric>
## ENSG00000003137.9    117.136       -2.31653  0.748442  -3.09514 1.96718e-03
## ENSG00000004468.13   440.636       -2.67812  0.684571  -3.91211 9.14919e-05
## ENSG00000004660.15   356.994       -2.16715  0.679544  -3.18911 1.42709e-03
## ENSG00000004776.13   309.782        2.60928  0.648061   4.02629 5.66648e-05
## ENSG00000005059.16   175.267        3.34492  0.583748   5.73007 1.00389e-08
## ...                      ...            ...       ...       ...         ...
## ENSG00000288825.1  983.26670        2.54303  1.008233   2.52226 1.16603e-02
## ENSG00000288841.1   27.57058        7.60671  1.542698   4.93078 8.19002e-07
## ENSG00000289690.1   22.52912        2.49508  0.715714   3.48614 4.90038e-04
## ENSG00000289692.1   56.62950        3.98232  1.417735   2.80893 4.97061e-03
## ENSG00000289697.1    4.90609        5.11596  2.364511   2.16364 3.04917e-02
##                           padj
##                      <numeric>
## ENSG00000003137.9  3.28501e-02
## ENSG00000004468.13 4.67995e-03
## ENSG00000004660.15 2.67659e-02
## ENSG00000004776.13 3.35615e-03
## ENSG00000005059.16 8.13394e-06
## ...                        ...
## ENSG00000288825.1  0.089432521
## ENSG00000288841.1  0.000176486
## ENSG00000289690.1  0.014235577
## ENSG00000289692.1  0.055935956
## ENSG00000289697.1  0.150161378

# Número de registros
nrow(filteredResults)

## [1] 1302

saveRDS(filteredResults, "./filteredResults.RDS")
write.csv(filteredResults, "./filteredResults.CSV")

# @see: "https://rdrr.io/bioc/ReportingTools/man/makeDESeqDF.html"
dfResults <- as.data.frame(filteredResults)

rm(resNotNans)

Añadir “gen name” y “gen change”

Posiblemente, por un problema con la versión de SALMON (según Antonio), la anotación de los IDs de los genes no permite asociarlos con el gen al cual referencian. Si se modifica la parte diferente, el problema se soluciona.

Ahora vamos a obtener una tabla que servirá para asociar el ID del gen con el nombre y símbolo del gen.

library(EnsDb.Hsapiens.v86)

## Loading required package: ensembldb

## Warning: package 'ensembldb' was built under R version 4.2.1

## Loading required package: GenomicFeatures

## Warning: package 'GenomicFeatures' was built under R version 4.2.1

## Loading required package: AnnotationDbi

## Warning: package 'AnnotationDbi' was built under R version 4.2.1

## Loading required package: AnnotationFilter

## Warning: package 'AnnotationFilter' was built under R version 4.2.1

## 
## Attaching package: 'ensembldb'

## The following object is masked from 'package:stats':
## 
##     filter

keytypes(EnsDb.Hsapiens.v86)

##  [1] "ENTREZID"            "EXONID"              "GENEBIOTYPE"        
##  [4] "GENEID"              "GENENAME"            "PROTDOMID"          
##  [7] "PROTEINDOMAINID"     "PROTEINDOMAINSOURCE" "PROTEINID"          
## [10] "SEQNAME"             "SEQSTRAND"           "SYMBOL"             
## [13] "TXBIOTYPE"           "TXID"                "TXNAME"             
## [16] "UNIPROTID"

columns(EnsDb.Hsapiens.v86)

##  [1] "ENTREZID"            "EXONID"              "EXONIDX"            
##  [4] "EXONSEQEND"          "EXONSEQSTART"        "GENEBIOTYPE"        
##  [7] "GENEID"              "GENENAME"            "GENESEQEND"         
## [10] "GENESEQSTART"        "INTERPROACCESSION"   "ISCIRCULAR"         
## [13] "PROTDOMEND"          "PROTDOMSTART"        "PROTEINDOMAINID"    
## [16] "PROTEINDOMAINSOURCE" "PROTEINID"           "PROTEINSEQUENCE"    
## [19] "SEQCOORDSYSTEM"      "SEQLENGTH"           "SEQNAME"            
## [22] "SEQSTRAND"           "SYMBOL"              "TXBIOTYPE"          
## [25] "TXCDSSEQEND"         "TXCDSSEQSTART"       "TXID"               
## [28] "TXNAME"              "TXSEQEND"            "TXSEQSTART"         
## [31] "UNIPROTDB"           "UNIPROTID"           "UNIPROTMAPPINGTYPE"

cols <- c("SYMBOL", "GENENAME")

# Por algún motivo, el ID del gen tiene una coletilla tipo: ".01". La base de datos
# no lo reconoce, así que es necesario eliminarlo
lTxID <- row.names(res)
lSplittedTxID <- unlist(strsplit(lTxID, "\\."))[seq(1, 2 * length(lTxID), 2)]

id2gene <- select(EnsDb.Hsapiens.v86, keys=lSplittedTxID, columns=cols, keytype="GENEID")

rm(lTxID, lSplittedTxID)

Añado el nombre del gen a los resultados. Aquí se pueden contabilizar los transcritos que no están asociados a ningún gen por ahora.

# Creo una función, porque volveré a utilizarla posteriormente
ObtainGenNames <- function(lGenIDs, id2gene) {
  genName <-
    unlist(lapply(lGenIDs, function(genid) {
      splittedGenId <- strsplit(genid, "\\.")[1]
      splittedGenId <- unlist(splittedGenId)[1]
      name <- id2gene[id2gene$GENEID == splittedGenId,]$GENENAME
      if (length(name) != 0) {
        return(name)
      } else {
        return("Not Available")
      }
    }))
  return(genName)
}

genName <- ObtainGenNames(rownames(filteredResults), id2gene = id2gene)

# Número de transcritos que no estan asociados a ningún gen...
length(genName[genName == "Not Available"])

## [1] 40

dfResults$`gen name` <- genName

rm(filteredResults)

Para facilitar la comprensión de los resultados, se va a generar una nueva columna en el dataframe de los resultados para indicar si se reprime o se sobreexpresa el transcrito en concreto.

# También se guarda esta función porque se va utilizar más adelante
ObtainChange <- function(lLog2FC) {
  # Aquí guardo unas variables por si se quiere cambiar el texto de la variable "change".
  sOverexpressed <- "overexpressed"
  sRepressed <- "repressed"
  sNoChange <- "no change"
  
  # Generar una variable que indica si se sobreexpresa o se reprime el transcrito
  # Se genera en base al valor de la variable "log2FoldChange"
  lGenChange <- unlist(lapply(lLog2FC, function(FC) {
    if (abs(FC) > 0) {
      return(sOverexpressed)
    } else if (abs(FC) < 0) {
      return(sRepressed)
    } else {
      return(sNoChange)
    }
  }))
  
  return(lGenChange)
}

# Añadir la nueva variable al dataframe con los resultados
dfResults$`change` <- ObtainChange(dfResults$log2FoldChange)

Guardamos el dataframe generado con los resultados.

# Guardar como un objeto serializado y como un csv
saveRDS(dfResults, "./finalResults.RDS")
write.csv(dfResults, "./finalResults.CSV")

Contabilizar genes sobreexpresados y reprimidos

En este apartado se contabilizan el número de genes que se sobreexpresan o se reprimen de forma más significativa. Es decir, los que pasan el filtro del p-valor y el Fold Change.

# Contabilizar el número de genes sobreexpresados y reprimidos
# Sobreexpresados
nrow(dfResults[dfResults$log2FoldChange > 0,])

## [1] 971

# Reprimidos
nrow(dfResults[dfResults$log2FoldChange < 0,])

## [1] 331

Volcano plot

Mostramos en forma de gráfica los resultados.

# Representando la misma gráfica pero ggplot2
library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.2.2

library(ggrepel)

## Warning: package 'ggrepel' was built under R version 4.2.1

# Para mostrar todos los puntos del volcano plot, es necesario construir un
# dataframe con toda la información

# Genero un nuevo dataframe a partir del resultado de DESeq
dfToPlot <- as.data.frame(res)
# Añadir el nombre del gen
dfToPlot$`gen name` <- ObtainGenNames(rownames(res), id2gene = id2gene)
# Añadir el "change"
lChange <- unlist(lapply(1:nrow(dfToPlot), function(iRow) {
  FC <- dfToPlot[iRow,]$log2FoldChange
  pv <- dfToPlot[iRow,]$pvalue
  
  if (is.na(FC) || is.na(pv)) {
    return("not available")
  }
  
  if (abs(FC) > log2FC && pv < pValue) {
    return("log2FC and p-value")
  } else if (abs(FC) > log2FC && pv > pValue) {
    return("log2FC")
  } else if (abs(FC) < log2FC && pv < pValue) {
    return("p-value")
  } else {
    return("no change")
  }
}))
dfToPlot$change <- lChange

# Parámetros de la gráfica
dLineWidth = 0.5

# Graficar
volcano <- ggplot(dfToPlot, aes(x = log2FoldChange, y = -log10(pvalue), color = change, text = paste0("gene: ", `gen name`, " codigo: ", rownames(dfToPlot)))) + 
  geom_point() + 
  scale_color_manual(
    breaks = c("not available", "no change", "log2FC", "p-value", "log2FC and p-value"), 
    values=c("dark gray", "grey", "#F8766D", "#619CFF", "#E76BF3")) +
  geom_hline(yintercept=-log10(pValue), linetype="dashed", color = "black", linewidth = dLineWidth) +
  geom_vline(xintercept=-log2FC, linetype="dashed", color = "black", linewidth = dLineWidth) +
  geom_vline(xintercept=log2FC, linetype="dashed", color = "black", linewidth = dLineWidth) +
  ggtitle("Volcano plot")

# Mostramos el gráfico con plotly para poder ver cada punto con que gen corresponde
plotly::ggplotly(volcano)

## Warning: `gather_()` was deprecated in tidyr 1.2.0.
## ℹ Please use `gather()` instead.
## ℹ The deprecated feature was likely used in the plotly package.
##   Please report the issue at <https://github.com/plotly/plotly.R/issues>.

rm(dfToPlot, volcano)